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TECHNICAL MEMORANDUM 


TRANSPORT PHENOMENA IN THE MICROPORES OF 
PLUG-TYPE PHASE SEPARATORS 

I. THE PROBLEM 


A. Introduction 

This study investigates the transport phenomena within a micropore of a porous-plug trapping 
device. The porous plug has been used as a phase separator in liquid helium systems and can potentially 
be used as a pressure control device as well. In ground- or space-based storage applications, heat 
invariably enters the tank causing the cryogen to boil off and evaporate, thus increasing the tank pres- 
sure. The tank pressure must be controlled to prevent tank over-pressurization and structural failure and 
to maintain the desired thermodynamic state for the particular application. Tank venting can be accom- 
plished easily in ground-based systems since the vapor location is predictable, namely at the top of the 
tank. In space-based systems where the acceleration is minimized, tank venting becomes more complex 
since the vapor location is unknown. Numerous studies in the past have examined both active and pas- 
sive techniques for venting and tank pressure reduction. 1 Active methods require complex systems that 
employ mechanical pumps and valves to achieve efficient venting of liquid-free vapor. Passive venting 
is much simpler since operation of the pressure control device is driven entirely by the pressure and 
temperature differential across the liquid/vapor interface and does not require the use of complex 
hardware. 

In the past, performance of porous-plug trapping devices has been estimated by assuming that the 
liquid in the individual pore is static and that the heat transfer through the pore is entirely by conduction. 
Recent work using pore sizes in the range of the porous plug, 0.1 to 1 .0 micrometers (/mi), have 
indicated that thermocapillary- and evaporation-driven convection may also be present. Therefore, the 
objective of this study is to investigate some of the physical parameters that influence the transport 
phenomena within the porous plug. 


B. Background 

Advanced cryogenic tankage requires the use of sophisticated insulation systems to minimize the 
heat into the fluid from the environment. Since the insulation cannot completely intercept the heat enter- 
ing the tank, a means of rejecting this heat is required. Heat rejection through venting is one common 
approach. 

Microgravity venting techniques studied in the past include both active and passive systems. 
Active systems offer more control over the venting process since the system can be designed to accom- 
modate variable conditions and high heat rejection rates. One technique for active venting in micro - 
gravity is to apply either continuous or intermittent thrust to settle the liquid into a predetermined loca- 
tion within the tank, thus allowing venting of liquid-free vapor. This operation complicates mission 
planning by requiring vent operations to coincide with RCS thrusting operations. Settled venting has 
been used as the primary technique on the Saturn S-IVB stage and is currently employed on the Centaur 
upper stages. Another active venting method that has been studied extensively is the thermodynamic 



vent system (TVS). 2 3 The TVS provides the combined function of heat extraction and mixing of the 
liquid and vapor, i.e., thermal destratification. Heat is rejected by passing a small amount of liquid 
through a Joule-Thompson (J-T) valve and a heat exchanger. When the liquid is passed through the J-T 
valve, it is flashed into a two-phase mixture, which results in a temperature decrease due to the pressure 
drop across the valve. This cooler fluid is then used as a heat sink by passing it through a heat exchanger 
which extracts energy from the circulating bulk liquid. The two-phase flow is superheated in the heat 
exchanger and subsequently vented as a gas. TVS concepts are highly efficient and can be designed to 
reject large amounts of heat, but they are relatively complex, heavy, and expensive. 

Passive systems based on porous-plug phase separation are much simpler and rely on a pressure 
differential and the thermodynamic properties of the cryogen to reject energy from the system. Porous 
plugs have been studied and tested previously with superfluid helium. 4 - 6 These devices mainly consist 
of a sintered metal or fine screen mesh and are difficult to analyze due to irregularly shaped flow pas- 
sages and inherent nonuniformities in the manufacturing process. These devices exploit the fluid's sur- 
face tension to form a phase boundary between the liquid and the vapor in the microscopic passages 

within the sintered metal or screens. These passages typically range from 0.1 to 1 .0 jim in diameter. 
Evaporation takes place across the liquid/vapor boundary to reject heat. An advanced porous plug that 
has been designed and tested recently 7 uses an advanced manufacturing process to fabricate a plug with 
multiple uniform cylindrical passages running parallel to its axis of symmetry. This type of design is 
easier to analyze due to its simple geometry. Figure 1 shows a typical wetted porous plug that is 
mounted in the wall of a cryogenic dewar or tank. 


Plug Material Pore Dia> 



■ Tank Wall 



Figure]. Porous plug. 

During nominal operations the plug provides a means of venting the tank and controlling tank 
pressure without expelling liquid across the pore surface. One means of passively controlling tank pres- 
sure and rejecting heat from a liquid is through evaporation across the liquid-vapor interface in a porous 
plug. If the liquid vapor pressure is greater than the ullage pressure, the liquid is superheated with 
respect to its saturated state. In order for the liquid to reach saturation, the liquid must reject heat or 
energy. This rejection is accomplished through evaporation within the porous plug. The porous plug can 
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perform this function regardless of whether liquid or vapor is in contact with the internal surface of the 
plug Past studies have shown that for the sintered porous plugs tested, there are three distinct observable 
flow regimes observed. 4 In the first regime, the pore is fully wetted, i.e„ filled with liquid, evaporation 
occurs across the interface, and the liquid assumes an equilibrium position as shown in figure 2. 1 he 
liquid/vapor interface residing on the vent side of the plug provides an evaporation site that transfers 
energy from the tank fluid and enables pressure control. The geometry of this interface is dictated by the 
jump momentum balance and includes the effects of pressure-, temperature-, and evaporation-induced 
recoil. In the second regime, the pore is dry and the vapor flow through the pore is laminar. In the third 
regime, the plug is also dry, as was the second regime, but the vapor flow is turbulent. This study con- 
centrates on the wetted case, namely regime one. 



Figure 2. Representative pore. 

An important aspect of the pore is the ability to sustain a pressure differential across the inter- 
face As shown in figure 1 , the pressure differential is caused by the tank pressure being greater than the 
vapor pressure on the exterior of the pore. The maximum pressure that the pore can sustain across the 
interface is commonly called the “bubble point” pressure and is theoretically described by the LaPlace 

Young equation: 


where D is the pore diameter. As long as the pressure difference across the plug is equal to or less than 
the value given by equation (1), the liquid will remain in the plug. Unlike sintered metal constructed 
porous plugs where the effective pore diameter is not constant and little control over the pore size is 
available, a porous plug with uniform cylindrical pores can be designed for a particular fluid and desired 

bubble point. 


C. Statement of the Problem 


In the wetted case (regime one) the heat rejection mechanism is much more complicated than in 
regimes two or three since the fluid dynamics in the liquid are more difficult to analyze than the vapor 
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flow through a pore. Heat must be rejected by evaporation across the interface. In the liquid, convection 
is driven by sidewall heating, evaporation, and thermocapillary effects. Most porous plug studies have 
assumed superfluid helium as the working fluid and only recently has hydrogen been discussed. 

In the previous studies, fluid dynamics within the pore could not be analyzed due to the compli- 
cated nature of the flow passages. Assumptions and simplifications were made to allow some under- 
standing of the fundamentals involved. The uniform pore design allows the use of computational fluid 
dynamics (CFD) techniques and can provide a more detailed investigation of the physical parameters 
that affect the transport phenomena within the pore. 

The objective of this study is to investigate the transport phenomena within the pores of a 
porous-plug trapping device. The effect of temperature differential across the interface and sidewall 
heating conditions has been examined. 


D. Approach 

The approach used to investigate the transport phenomena in a pore was to use CFD to model a 
representative pore. A reference finite element model developed to analyze a liquid acquisition device 
(LAD) pore , 8 described in section III, has been modified to account for the physical differences of the 
LAD and porous plug. The modified model is then used to investigate the influence of various physical 
parameters on the flow field within the pore and the mass transfer or evaporation across the pore inter- 
face. 


E. Summary 

Pressure control for space-based systems becomes more complicated than ground-based due to 
the uncertainty of the ullage location. Pressure control devices that achieve liquid-free venting have been 
studied in the past and include the porous-plug phase separator. Past studies had identified three flow 
regimes within the porous plug through experimentation. The complicated geometry and flow passages 
of the porous plugs have precluded detailed analysis until uniform cylindrical microscopic flow passages 
were able to be manufactured through advanced processes. This study numerically investigates the 
transport within and across a single cylindrical pore of porous plug that is wetted with liquid hydrogen. 


II. REVIEW OF THE LITERATURE 


A. Introduction 

The concept of phase separation using a porous media or plug has been studied in the past by 
several researchers. Past research is discussed in this section, namely, work that has been done using 
superfluid helium or He II, normal fluid helium, and initial work that has been done for liquid hydrogen 
using a porous plug as the phase separation device. Other work that pertains to the phase separation 
across an LAD pore is also discussed. The porous plug and LAD applications are different, but the 
governing equations for the flow within and across the LAD are also valid for the porous-plug flow 
field. 
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B. Porous-Plug Research 


As mentioned earlier, most of the past studies of sintered metal porous plug operation have been 
based on He II. Cooling requirements for experimental electronics have driven the need for He II phase 
separation devices. Karr and Urban 4 and Frederking et al. 5 investigated using the porous plug to modu- 
late helium vapor flow for the cooling of an experiment. Karr and Urban used an AI2-O3 ceramic plug 

with a nominal 5-/zm diameter pore size to evaluate both the steady and unsteady porous plug operation 
with helium. Steady-state data were obtained by measuring the mass flow rate through the plug as a 
function of the downstream pressure, and it was observed that the relation of m to AP shows two dis- 
tinct regimes, characterized by a change in slope of the curve. To experimentally evaluate the unsteady 
operation of the plug, a heater was installed downstream of the plug and used to modulate the vapor flow 
across the plug. This technique was effective in supplying cold helium vapor on demand. Frederking et 
al. used a flow control plate and shutter assembly to increase or decrease the evaporative area of the 
porous plug by opening or closing the shutter. This technique was also effective in regulating the supply 
of cold vapor. In both studies the standard two-fluid model for He II is used to theoretically characterize 
the mass flow through the pore. 

Hendricks and Karr, 6 7 9 have investigated the flow regimes previously identified in reference 4 
and the transition points between the regimes. In reference 6, a 316 stainless steel plug with an 0.5- /tm 
filtration grade was used to experimentally examine the transition between the lower and upper regimes. 
The lower regime was proposed to be defined when the liquid/vapor interface is downstream of the plug 
and the plug is filled with liquid. It was shown that the regime, transition occurs at the critical mass flow 
rate, m crit . In the upper regime, it was proposed that the liquid/vapor interface retreats into the plug, thus 
filling the pores with vapor. One explanation for this phenomena is that as the heat flux on the plug is 
increased, there comes a point where the liquid can no longer transfer the heat to the interface; therefore, 
the interface must move such that the energy is balanced. 

Schotte 10 examined He II phase separation using a narrow annular slit instead of the sinter metal 
porous plug. As with references 6, 7, and 9, Schotte measured the relation of m to APand observed the 
same transition occurring at the critical mass flow rate. Using the two fluid model, Schotte was able to 
show a good comparison with the experimental data. This was probably due to the uniformity and sim- 
plicity of the experiment geometry. 

Recently, utilization of a porous plug with fluids other than He II have been evaluated. Hendricks 
and Dingus 1 1 used normal fluid He or 4 He to evaluate the bubble point of a porous plug constructed of 
fc 2 -0 2 or jewelers rouge, as a function of the 4 He liquid temperature. The test data showed that the 
observed bubble point for 4 He decreased with increasing liquid temperature; therefore, the liquid must 
be precooled prior to entering the plug to prevent retention loss. Of more direct interest is work done by 
Hendricks and Dingus 12 on using a porous plug in a liquid hydrogen system. The porous plug used con- 
sisted of multiple uniform pores that were parallel with the plugs axis of symmetry. The advantage is 
that it allowed for more detailed analysis of the flow field within the individual pore since the geometry 
can be modeled more accurately. To evaluate the evaporation across the interface of the pore, a simple 
one-dimensional fin model was used. The analysis assumed that heat is transferred to the interface by 
solid conduction through the plug material and that liquid convection within the pore interior is negli- 
gible. 
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C. Relevant Numerical Models 


Sasmal and Hochstein 13 developed a finite volume model to predict the transient behavior of the 
thermocapillary convection within a two-dimensional cavity with a deforming interface. The fluid 
motion, thermocapillary convection, is driven by a variation in surface tension due to an imposed tem- 
perature gradient. In this case, the left boundary was held at a higher temperature than the right boundary 
and zero heat transfer was specified across the free surface and the lower boundary. Since the condition 
of zero heat transfer across the surface was imposed, no evaporation or interfacial mass transfer was 
considered. Sasmal and Hochstein modeled various free-surface geometric configurations including flat 
surface; acute contact angle, i.e., wetting fluid; and obtuse contact angle, i.e., nonwetting fluid. The 
obtuse contact angle case is more representative of the problem considered in this study since it repre- 
sents the surface due to the pressure differential between the tank ullage and the exterior of the plug. 
Their analysis showed that the free surface rises along the cold wall and falls along the hot wall. This 
phenomena was most likely due to the fluid being driven away from the hot and toward the cold bound- 
ary by the thermocapillary stress condition. 

Schmidt 8 developed a finite element model to investigate the steady-state behavior of the flow 
about an evaporating or condensing LAD pore, while accounting for thermocapillary convection. The 
problem domain consisted of a two-dimensional cavity having a curved free surface with an acute con- 
tact angle. A temperature gradient was imposed across the free surface to examine the influence of inter- 
facial mass transfer and thermocapillary convection on the liquid retention of the pore. This model is 
discussed in more detail in section III. One of the more significant findings by Schmidt was that thermo- 
capillary convection arising from condensation can cause a retention loss, while thermocapillary con- 
vection arising from evaporation causes a decrease in the surface curvature and, therefore, tends to make 
the interface remain stable. 


D. Summary 

In this section, previous experimental and analytical work pertaining to the porous plug trapping 
devices is discussed. Studies which concentrated on He II, normal He, and preliminary work with liquid 
hydrogen was summarized. Also, numerical models that are relevant to the problem addressed in this 
study are discussed. 


III. DESCRIPTION OF REFERENCE MODEL 


A. Introduction 

As shown in section II, Schmidt 8 has developed a two-dimensional finite element model of 
evaporation or condensation within an LAD pore. The model accounts for thermocapillary convection at 
the pore interface and the associated convection within the pore along with the mass transfer at the inter- 
face. The pore’s free surface is modeled with respect to the physical parameters at the interface, such as 
temperature, pressure, and fluid velocity. This model is used as the basis for the analysis of this study. 
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B. Physical and Numerical Model 

The physical domain of the reference model approximates that of the cross section of a single 
two dimensional pore of the LAD. The cross section, shown in figure 3, shows an individual interface or 
meniscus attached to individual circular screen wires. 

Interface Vapor 



Figure 3. LAD pore schematic. 

The boundary conditions for this representation are complex. The curved surface representing the 
screen wires and the free flow or unspecified flow boundaries at the left and right side of the domain 
would make this problem complicated to solve and, therefore, would require a much more sophisticated 
numerical model. In order to simplify the problem, several assumptions were made to obtain the prob- 
lem domain. The first assumption is that the interface diameter is small compared to the screen wire 
diameter; therefore, the curved screen boundaries can be approximated as flat. A second simplifying 
assumption is that the sidewall boundary conditions are no-slip. Figure 4 shows the problem domain that 
is used in the reference model. 



xl =0 xl=D 


Figure 4. Reference model problem domain. 
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The governing equations for the reference model are based on constant properties and derived 
from the continuity, momentum, and energy equations for an incompressible fluid, as shown below. 


Vij = o 


( 2 ) 


dV: 


p ~dt + pVjViJ + Pd ' [ ' -pP( t ~ T > h = 0 . 


(3) 


where 


and 


— + v T - CT = 0 

J j ?,jj u ’ 


(4) 



(5) 




(6) 


The Boussinesq approximation is used in the momentum equation (3) to account for buoyant 
forces. This also restricts the thermal dependence of density as a body force. 

A fourth governing equation defines the free surface shape due to the flow and temperature field 
characteristics. The normal component of the jump momentum balance on the surface yields a relation 
that represents the surface curvature and is given by: 


■2 


1 * = - I’d - r h ) + j(f p - 1 ) + 2 WjHjn, , 


(7) 


where k is the surface curvature, P {v) is the vapor pressure, P d is the liquid dynamic pressure, P h is the 
hydrostatic pressure,; is the evaporation flux, f p is the ratio of the liquid density to vapor density, and 
term 2/jVj jn J n i represents the viscous stress on the surface. 

The surface curvature can also be defined using the expression for the radius of curvature in two 
dimensions. 


K = - 


1 d 

dx 



( 8 ) 


Two additional relations are used, one to account for thermocapillarity and the other to account 
for the mass transfer across the interface. In order to account for the thermocapillary effect of the fluid, 
i.e., the thermal dependence of surface tension, a linear relation of surface tension as a function of tem- 
perature is used: 


7 = Ti- 


dy 

~dr 


\(T~ T ,). 


(9) 
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To account for the evaporation across the interface, the relation given below 14 is used to relate 
the mass transfer to the liquid temperature and its thermophysical properties: 


f 

j=<* 

\ 


2/tA/yv j 


1/2 


ip-Po) . 


( 10 ) 


where a is the coefficient of evaporation that estimates the fraction of impinging liquid molecules that 
do not evaporate, i.e., some of the molecules rebound off the surface rather than evaporating. The maxi- 
mum evaporation flux occurs when a = 1 . The universal gas constant and the molar weight of the vapor 
are given by R g and M w , respectively. The density difference in equation (10) represents the divergence 
from saturation to a superheated state corresponding to temperature T. It was shown 8 that the mass flux 
can be related to a temperature difference across the interface and takes the form of: 


J = 


a po V) L 

■ 0/2 


M,. 


\>'2 


\ 2nR s ) 


( [T-T 0 ) 


(ID 


All the aforementioned equations and relations are then nondimensionalized using the following 
scaling relations: 


x, = Dx, 


V ; = — V* , 
' D 


t = 



gi= a iSo . 


pv 2 


O:; = 

lJ D 


2 °ij 



( 12 ) 


j = 


pv .* 

T J ' 


T = T*\AT\ + T 0 . 


Using the scaling relations given in equation (12), the following nondimensional governing 
equations are obtained: 


where 


Vu= 0 , 


dv [ 

dt 


+ VjVlj + P*d,i-'ijJ- Gra i T *= 0 ’ 


dT * * 1 * 

— + V :T : - — T :; = 0 , 
dt J ' J Pr ' JJ 


Gr = 


v 2 


Grashoff number 


(13) 

(14) 

(15) 

(16) 
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Prandtl number 


(17) 



Using equation (7), the equation for the interfacial curvature is given by: 


* 

K = 


_ fl 0 y W * - Ca(p* d - P (v) * ) + V T* 2 + 2CaV* j n j ,«,• 


where. 


and from equation (8), 


1 - CrT 


Bo = 


Ca - 


PgDj 

Yi 

pv 2 

Y/D 


Rs 

dy 


v r= 2 

Rs 2 


Cr — — 

Y 


dr 


|A71 , 


bond number 
capillary number 

recoil parameter 
crispation number 


k =-- 


J d_ 

y\f dx |L 


r r y ,21 

-1/2] 


i 


Using equation (11), the relation for the mass transfer across the interface becomes: 


.* T 
J ~ Rs ’ 


where 


(18) 

(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 


/ p vr o 3/2 [2 nR g 


Rs 


1/2 


_ J P 


aDL\AT\ M, 


evaporation resistance . 


w 


(25) 


Boundary conditions must be imposed on the four boundaries shown in figure 4, the right and 
left sides and the upper and lower surfaces. Making the assumption of no-slip and impermeability on 
both the left and right sides, and making the further assumption that the sidewall temperature is specified 
and constant, the conditions on boundaries 1 and 3 are defined, namely: 



( 26 ) 
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( 27 ) 


V 2 *=0 , 
and 

T* = 7’, . (28) 


The boundary conditions on the upper and lower surfaces are more complex. The lower surface 
must satisfy the restriction that the flow entering or exiting the cavity through the lower surface be uni- 
form and equal to the flow leaving or entering the cavity through the upper surface due to evaporation or 
condensation. Also, the temperature on the lower surface is assumed to be constant and equal to the side 
boundary temperature. The lower boundary conditions are: 


V, =0 , 


(29) 


and 


v* 2 =— J ^T 4 - 

2 r 2 Jr 4 Rs 


T =T, . 


(30) 

(31) 


For the upper free surface, one must consider the normal and tangential stress conditions due to 
the surface curvature and the thermocapillarity, respectively, and, in addition, one must account for 
liquid ejected from or deposited onto the surface due to evaporation or condensation. The upper bound- 
ary conditions are: 


* ^ _* A f(l 

T unj = -ReTj=-— Ti , 


v 2 = — 


2 

n 2 


( r r* ^ 

k R S , 


‘1 


(32) 

(33) 


where 


and 


r ( n,- = - NuT 


(34) 


Re = 


dy 


dr 


D|A7j 


Ma = 


pv 2 

dy 


dr 


£>|A7j 


Nu = 


hD 

t(0 


surface tension Reynolds number 


Marangoni number 


Nusselt number . 


(35) 


(36) 

(37) 


An alternate definition for the Nusselt number, which gives the rate of latent heat transfer relative 
to the heat conduction is termed the mass transfer Nusselt number and is given in equation (38). It 
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should be noted that in reference 8, equations (37) and (38) are referred to as the convective Biot num- 
ber, Bi c , and the mass transfer Biot number, Bi m , respectively. 


Nu = ■ 


RsE 


mass transfer Nusselt number 


(38) 


where the nondimensional group, E, shown in equation (38) that was not previously defined is given by: 




* 1*71 

pvL 


evaporation number . 


(39) 


As discussed in reference 8, both representations of the Nusselt number, equations (37) and (38), 
are equivalent in terms of defining the temperature distribution along the surface. 

The surface solution, y*> is determined by equating equations (18) and (23) and solving using the 
following boundary conditions: 


= 1 , at jc 7 = 0 and 1 , 


dy^_ 

dx* 


= tan(7r / 2 - 0 )) 


at xl = 0 and 1 . 


(40) 

(41) 


C. Summary 

The reference model described in this section solves for the flow field and thermal profile in a 
simplified liquid acquisition device pore. The model employs a finite element numerical method to solve 
the governing equations and boundary conditions for velocity, temperature, pressure, and surface posi- 
tion. The model accounts for thermocapillary and evaporation or condensation-driven convection within 
the cavity. 


IV. MODIFICATIONS TO REFERENCE MODEL 


A. Introduction 

The problem considered in this study, the porous plug, is similar in many ways to the problem 
solved using the reference model. 8 Many of the assumptions still apply, namely, incompressibility, con- 
stant properties, and limiting the domain to two dimensions. Therefore, the governing equations are 
valid, i.e., continuity, momentum, and energy. Some of the boundary conditions of the domain must be 
modified to accurately model the problem addressed in this study. 

The modifications that are required are the calculation of the free surface shape, the specification 
of the fully developed flow condition on the bottom boundary, and the ability to apply a temperature 
gradient on the side boundary as well as an isothermal boundary. The aspect ratio of the problem is also 
increased from one to two, but is considered minor and will not be discussed in detail. 
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B. Free Surface Boundary Condition 


As shown in section I, the free surface for a fully wetted porous plug differs from that of the 
LAD pore due to the direction of the pressure differential across the interface. In the LAD application 
the ullage pressure is greater than that of the liquid; therefore, the pressure differential tends to drive the 
ullage and the interface into the liquid. In the fully wetted porous plug application, the liquid pressure is 
greater than the ullage pressure; therefore, the reverse is true, namely, the pressure differential tends to 

drive the liquid and interface into the ullage. This fundamental difference in the direction of the Ap 
requires that the surface solution be modified. Although the surface solution must be modified, the tan- 
gential stress condition, i.e., the thermocapillary effect, is still valid. 

In order to define the modified surface, a change in the coordinate system is required. The coor- 
dinate system used in the free surface solution is shown in figure 5: 

x2, y 



xl=x= 0 xl=x= 1 

Figure 5. Free surface coordinate system. 

where v 1 = .v, v = v^-.v2 and v* v) is the y coordinate of the surface at each x. In order to calculate the free 
surface, equations ( 1 8) through (23) are used. Expanding the derivative in the right-hand term, equation 
(23) can be written as: 


* 

K 



(42) 


Equating equations (18) and (42) gives: 




1 + 


W 2 


- 3 / 2 " 


g„y (,) * + n 


1 - D 


(43) 
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where: 


n = -Ca(p* u +P* d - P {v) * ) + V r T* 2 , 
and 

D - \- CrT* . 


(44) 

(45) 


In equation (44), the viscous stress term, 2CaV* jitjitj, is dropped since its influence on the sur- 
face curvature is negligible 8 and the term P*, the tank ullage pressure, is added to account for the higher 
liquid pressure in the porous plug. 

Rearranging equation (43) and dropping the * notation, the following second-order, nonlinear 
differential equation is obtained: 


y,xx ~ 


f r> - (s ) _j_ j-j \ 


Boy K 


1 - D 


K4' 2 =° • 


(46) 


The boundary conditions on equation (46) consist of a Dirchlet and Neumann condition. The 
Dirichlet condition is that at x = 0 and x = 1 , y = 0. This is commonly referred to as a pinning condition. 
Unlike an LAD, the boundary conditions at the end point are independent of the fluid contact angle. The 
Neumann condition restrains the slope of the interface at the center line to zero. This also implies that 
the interface is symmetrical about the center line. 

In order to solve for the interface shape, a finite difference representation of equation (46) is 
obtained: 


v l+ i ~ 2y,- + y,-_ i 

h 2 


Boy i + n 
1-D 


! | ( y,-+i -y,-_i 
2h 


\2 

) 


-| 3/2 


= o 


(47) 


where h is equal to the separation distance of the nodes. 

Equation (47) is solved using Newton’s method for a nonlinear series 15 to converge on a solution 
for the independent variable y. Newton’s method uses the following iteration algorithm: 


where: 


y(k) = y (k-D _ A (k-l) 


AV 


(*-b 


('rT 1 ^- 0 * 


(48) 

(49) 


and where J tj is defined as the Jacobian of the root vector, Fj. 

Once both the Jacobian and the root equations are determined, the surface solution is found by 
iteration. Starting with an initial guess for y/, the value of is determined and a new value of yi is 

obtained using equation (48). This process is repeated until, yf k) - y^ k ~ ]) < stol, where stol is the surface 
convergence tolerance and is equal to 0.001 in this study. 
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Recalling equation (46), it is seen that the nondimensional terms Ca, V r , and Cr are defined by 
the fluid under consideration, and the values of dynamic pressure and temperature are determined from 
the solution of the fluid governing equations. The parameter that can be specified directly is 
z\p* — — p * ' ** . This parameter can be used to investigate the influence of the differential pressure 

across the pore. 

One limitation of this surface model is that as the AP * approaches A P* b , the bubble point of the 
pore, the solution is indeterminate. This limitation occurs as the slope of the surface at x = 0 and 1 
approaches infinity. Figure 6 shows the effect of AP on the free surface shape. 


2.40 
2.30 
2.20 
2.10 
2.00 

0.0 0.2 0.4 0.6 0.8 1.0 

xl 

Figure 6. Free surface shape versus pore pressure differential. 



C. Lower Boundary Condition 

As described in section III. the reference model was developed to analyze the parameters associ- 
ated with the loss of retention in a LAD. In that problem, the lower boundary condition was assumed to 
be one of uniform flow across the boundary, with the total flow equal to the amount evaporated or con- 
densed. For the case of evaporation, this assumes that the fluid leaving the cavity is being replaced from 
a large reservoir of constant temperature beneath the interface. For the LAD, a uniform flow boundary 
condition is valid, however, this boundary condition would not be valid for the porous plug. Physically, 
a single pore in a porous plug has a large aspect ratio, i.e., the ratio of its length to its diameter, on the 
order of several hundred. In this study, a pore length of two diameters is used in order to reduce the 
computation requirements. Therefore, a lower boundary condition which represents a fully developed 
flow is more appropriate. 

To develop the lower velocity boundary condition, the expression for fully developed laminar 
flow in a channel 16 is used and is shown below: 


u~ u 


max 



(50) 
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Equation (50) defines a parabolic velocity distribution u(x) that has a value of zero at the 
sidewalls, x = ±a, and is maximum at the center line, * = 0. The variable « max is the maximum value of 
the velocity distribution, a is the half channel width, and x is the independent variable. Assuming the 
domain shown in figure 7, equation (50) can be written as shown in equation (51). 


y 


\ 


X 

\* — D-2a — 



Figure 7. Lower boundary condition coordinate frame. 

(x-a ) 2 


u(x ) = w n 


a 


(51) 


Using the scaling relations defined in equation (12), equation (51) is nondimensionalized and is 
given by: 


* 

u 


* 

— ^max 



* 

X 



(52) 


In order to find the total mass flow, J* , across the boundary, equation (52) is integrated from 0 
to 1 . ‘ 


J* = ju*dx* = jw’ 




0 0 

and solving for w max the following expression is obtained: 


<fr*=!«max 


where 


M max 2 ^ 


J* = J/ dr 4 


(53) 


(54) 

(55) 
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Using equations (52) and (54) along with the total mass flux that is exiting the pore through 
evaporation, the lower boundary condition is defined. A typical velocity profile is shown in figure 8 



Figure 8. Lower boundary velocity distribution. 


D. Sidewall Boundary Condition 

In the reference model, the sidewall boundary conditions were assumed to be isothermal. The 
isothermal boundary condition was appropriate for the problem studied 8 because of the small length 
scale between the liquid and the vapor; however, for the porous plug, the condition is valid only if the 
porous plug material is manufactured from a high thermal conductivity material such as copper. If the 
porous plug is manufactured from a lower thermal conductivity material, the isothermal boundary con- 
ditions would not apply, rather there will be a thermal gradient through the material. For this case, a 
linear temperature gradient will be imposed on the sidewall boundaries of the porous plug, where the 
magnitude will be based on the plug thermal conductivity. 

Figure 9 shows a linear temperature profile along the side boundaries of the porous plug. 



Figure 9. Side boundary temperature distribution. 
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In order to model the thermal gradient boundary condition, an assumption that the temperature 
distribution is linear and can be approximated as a one-dimensional conduction problem is made. 
Assuming that the temperature distribution through the solid boundary is linear implies that the thermal 
conductivity of the material is constant, i.e., the thermal conductivity is independent of temperature. An 
additional assumption that there is no heat generation from within the plug material is also made. With 
the preceding assumptions, the equation for the temperature through the plug 17 is given by: 


5S 

II 

o 

@ 0 < y< L 

(56) 

with both Dirchlet and Neumann boundary conditions given by: 



^ = T'liq * 
and 

o 

II 

(57) 

kT y + h(T-T vap ) = 0 . 

II 

® 

(58) 

Using the appropriate scaling factors given in equation (12), equations (56) through (58) can be 
nondimensionalized as: 

o 

II 

@ 0 < y* < 1 

(59) 

T* = 1 , 

O 

II 

* 

@ 

(60) 

T* y + BiT* = 0 , 

@y* = AR 

(61) 


where the aspect ratio, AR, is defined as the ratio of the numerical domain length to the pore diameter 
and where the Biot number, Bi, is defined as the ratio of the heat transfer by convection to that by con- 
duction across a material of length L. A small value of Bi implies that the heat transfer is dominated by 
conduction. Based on this definition, the Biot number takes on the form: 


Bi = 


hL 

fcisol) ’ 


(62) 


where h is the convective heat transfer coefficient based on the vapor, k is the thermal conductivity of 
the plug material, and L is the physical length of the plug. 

Integrating equation (59) twice and applying the boundary conditions shown in equations (60) 
and (61) gives the following relation for the sidewall temperature distribution: 


T = 1 - 


Bi 


1 + (AR)Bi 


-x2 


(63) 


Equation (63) is used as the left and right wall temperature boundary conditions for the case of 
nonisothermal sidewalls. Figure 10 shows the effect of Bi on the temperature distribution through the 
plug material. 
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Figure 10. Sidewall temperature for various Bi. 


E. Summary 

In this section, the modifications to the existing model that were required to analyze the porous 
plug are discussed. The only modifications that were required included some of the boundary conditions, 
namely, the free surface solution, the redefinition of the bottom velocity boundary condition, and the 
right and left wall temperature boundary conditions. 


V. RESULTS AND FINDINGS 


A. Introduction 

Using the model discussed in section IV, the porous plug temperature profile, flow field, and 
mass transfer across the interface are determined for various initial conditions. In this section, the 
hydrogen properties and associated nondimensional parameters that are used in this study are presented 
and discussed. Also, a parametric analysis on the model grid size is discussed for the case of combined 
evaporation/thermocapillary flow. This is done to evaluate the effect of grid refinement on the numerical 
solution and to choose an appropriate grid resolution for the analysis. 

Using the hydrogen parameters identified and the appropriate grid size, three cases are discussed. 
(1) the basic state, pure evaporation driven flow; (2) combined evaporation/thermocapillary flow; and 
(3) combined evaporation/thermocapillary flow with a sidewall temperature gradient. Cases (2) and (3) 
are run for various A T magnitudes across the pore and are compared to case (1 )• 


B. Fluid Properties and Nondimensional Groups 

As was discussed in section II, most prior work in the area of porous plug research has been 
based on helium as the working fluid. Only cursory examinations have been done using liquid hydrogen. 
In this section, the relevant hydrogen properties that are used in this study are given in table 1 , and the 
associated nondimensional groupings are calculated and presented. The hydrogen properties are taken at 
ambient pressure and saturated conditions. 


19 




Table 1 . Saturated hydrogen properties. 



Units 

Value 

P sat 

(MPa) 

0.1013 

T sat 

(K) 

20.3 

p(D 


70.79 

p {x) 

(1) 

(kg/m 3 ) 

1.339 

tl K) 

(kg/m-s) 

1.319X10- 5 

n {v) 

(kg/m-s) 

8.5x1 0- 7 

y(0 

(m 2 /s) 

1. 860xl0- 7 

V (v) 

(m 2 /s) 

6.348xl0- 7 

k U) 

(W/m-K) 

1.182x10-' 

k (v) 

(W/m-K) 

1.883X10- 2 

Cp {l) 

(J/kg-K) 

9.689X10 3 

Cp (v) 

(J/kg-K) 

1.176X10 4 

a (t) 

(m 2 /s) 

1.440x10-7 

a (v) 

(m 2 /s) 

1.042x10-* 

L 

(J/kg) 

4.456x10 s 

M w 

(kg/kgmol) 

2.016 

Y 

(N/m) 

1.973x10-3 

\dy/dT\ 

(N/m-K) 

1.620X10- 4 

P 

(1/K) 

1.569X10- 2 


(J/kmol-K) 

8.314X10 3 


Using the property values given in table 1 , the corresponding relevant nondimensional groups 
that are independent of the temperature differential, AT, across the pore are given in table 2. For all non- 
dimensional parameters, a pore diameter of 1 /rm and a coefficient of evaporation of 0.5 are used. 

Table 2. Nondimensional groups. 


Parameter 

Value 

Nu 

7.639x10' 

Bo 

3.52x10-7 

Ca 

1. 241xl0- 3 

fp 

5.287x10' 


1.911x10-' 

A 

1.288X10 7 

Pr 

1.292 
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Based on the value of Ca given in table 2, the nondimensional bubble point, A P b = UCa , is 
equal to 1 ,6 1 1 . 

Table 3 gives the relevant nondimensional groups that are dependent on AT. Also, for the case 

where Ma is nonzero, Ca can be defined as CrPr/Ma and the nondimensional bubble point, A P p = 
2Ma/CrPr , is equal to 1,608. In this study, the effect of A P across the pore is not examined, but is set at 
a fixed value of 1 ,000. 


Table 3. Nondimensional groups versus AT. 


Parameter 

© 

II 

< 

Value 
AT = 1.25 

in 

ll 

<3 

Cr 

8.21 lx 10- 2 

1.026x10-' 

1.232x10-' 

E 

2.015x10-2 

2.518X10- 2 

3.022x10-2 

Gr 

4.449x 10- 6 

5.561XKH 

6.674X 10 I- 6 

Ma 

8.529x10' 

1.066X10 2 

1.279X10 2 

Rs 

6.498x10-’ 

5.199x10-’ 

4.332x10-' 

Vr 

9.908X10- 2 

1.238x 10-’ 

1.486x10-' 


C. Grid Refinement 

In order to determine the appropriate grid size to use in this study, the case of combined evapo- 
ration/thermocapillary flow is used. For this grid refinement exercise, the values of the nondimensional 
groups used are given in tables 2 and 3 for AT = 1 . For this case, four grid sizes were chosen where the 
number of elements in the x 1 direction is given by MXE , and the number of elements in the .\‘2 direction 

is given by NYE. Grid sizes of I Ox 15^25x30 in increments of five elements for NXE and NYE are used, 
where the total number of grid points is equal to (2* NXE+\)*(2*NYE+\). The grid size of 10x15 corre- 
sponded to the coarsest grid required for adequate velocity convergence, and the grid size of 25x30 cor- 
responded to the finest grid that could be solved using the available computer resources. 

To evaluate the effect of grid size, the surface temperature distributions, the mean surface tem- 
perature, and the evaporation flux for the four grid sizes are compared. The surface temperature is used 
since it drives the entire flow field, either by causing evaporation across the interface or by causing 
thermocapillary convection on the surface. 

For the case of combined flow with a AT = 1, the surface temperature gradient for various grid 
resolutions is given in figure 1 1 . 

The magnitude of the gradient and the percent change due to grid refinement is shown in table 4. 

In addition to the surface temperature distribution, the evaporation flux is also used to evaluate 
the grid size. Table 5 gives the percent change in J* due to grid refinement. 
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Figure 1 1 . Surface temperature for various grid sizes, AT= 1 . 
Table 4. Surface temperature gradient sensitivity to grid size, A T= 1. 


Grid Size 

Tjc 

Percent 

Change 

10x15 

-17.98 

- 

15x20 

-25.67 

42.8 

20x25 

-32.62 

27.1 

25x30 

-39.05 

19.7 


Table 5. Evaporation flux sensitivity to grid size, AT = 1. 


Grid Size 

J* 

Percent 

Change 

10x15 

0.077 

- 

15x20 

0.082 

6.5 

20x25 

0.086 

4.9 

25x30 

0.089 

3.5 


As is shown in both table 4 and 5, the increased number of elements or grid points serves to help 
resolve the large temperature gradient in the pore comers and in the overall evaporation flux across the 
pore. 


Although a grid size of 25x30 gives approximately a 20-percent change in T* from the previous 
grid, the value of J* changes by only 3.5 percent. Therefore, due to the relatively small change in J* and 
due to the exorbitant computational requirements of the 25x30 grid, a grid size of 20x25 is used in this 
study. It should be noted that a more refined grid should help to resolve the comer temperature gradient 
further and produced more accurate results in J*. 
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D. Basic State (Pure Evaporation) 


For the purpose of this study, the basic state for the porous plug analysis will be defined as pure 
evaporation across the interface, i.e., setting the Marangoni number equal to zero with a AT = 1.0 K. All 
other nondimensional parameters are as shown in tables 2 and 3 for A T = 1.0 K are used. Also, the side- 
wall temperature boundary condition is taken as isothermal having a value of 1 .0. 

Figure 1 2 shows the temperature and stream function profile that results due to the latent heat 
transport across the surface due to evaporation. The stream function profile indicates that the liquid is 
being driven toward the comers of the pore due to the increased evaporation at the regions of higher 
surface temperature. 



Figure 12. Basic state temperature and stream function distribution. 

Figure 13 shows the surface temperature distribution for the basic state. As shown, there is a 
large surface temperature gradient in the comer regions of the pore. 

Based on this surface temperature distribution, the mean surface temperature, T s ur j , is equal to 

0.088, or in dimensional form 19.38 K, and total mass transfer or evaporation, J , across the surface is 
calculated using equation (24) and is found to be equal to 0.086, or in dimensional form, 1 .1 36 kg/m 2 -s, 
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Figure 13. Basic state surface temperature distribution. 


E. Combined Evaporation/Thermocapillary Flow 

The case of combined flow is indicative of an actual pore due to the presence of thermocapillary 
convection. This section assumes combined flow while specifying isothermal sidewall boundary 
conditions, i.e., Twaii = 1.0. As can be seen in figure 14, no significant change in the isotherms is 
observed, but considerable change in the flow field is apparent. The flow field for combined flow is 
characterized by two cells of equal and opposite stream function distributions. This flow field results 
from the liquid on the surface being driven into the center of the pore due to the variation of surface 
tension. Figure 1 4 shows the temperature and stream function distribution for combined flow with an 
imposed A T = 1 .0 K. 



x2 1 f 4 


T -1 


_J 1 I I t I I 


0 .2 
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Figure 14. Temperature and stream function distribution, Ma = 85.29. 
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For increasing levels of superheat, i.e., A T, Marangoni number will also increase due to its 
dependence on AT, thus increasing the liquid convection from the hot boundary, namely the sidewalls, to 
the center of the cavity. This liquid convection from the sidewall to the center will tend to increase T sur j 
by transporting warmer liquid along the surface toward the center of the pore. However, this increase in 
surface temperature will result in an increase in the evaporation across the interface, which tends to 
decrease T surf due to the latent heat transport from the liquid to the vapor. In addition, the increased 
evaporation will tend to drive the liquid to the comers, as evidenced in figure 12. Table 6 gives the 
values for T sur j and J* for various levels of super heat. 

Table 6. T sur j and J* versus A T. 


AT (K) 

T 

1 surf 

J* 

Basic 

0.088 

0.086 

1.0 

0.087 

0.086 

1.25 

0.086 

0.107 

1.5 

0.086 

0.129 


Using the scaling relations given in equation (12) and the relevant properties given in table 1 , the 
dimensional values of Tsurf and J are given in table 7. 

Table 7. T sur j and J* versus A T (dimensional). 


AT (K) 

Tsurf (K) 

J (kg/m 2 -s) 

Basic 

19.38 

1.136 

1.0 

19.38 

1.132 

1.25 

19.16 

1.410 

1.5 

18.93 

1.699 


Figure 15 show the variation of the pore mean surface temperature, T sur j, and evaporation flux, 
J *, respectively. 
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Figure 15. T surf and J* versus A T. 
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Although T sur f is nonlinear over the range of AT considered, J* is linear since it depends on the 

actual surface temperature, 7*, and little variation of the actual surface temperature is observed. It was 
shown by Schmidt 8 that surface temperature is strongly dependent on the Nusselt number, Nu , which is 
independent of AT. It was also shown that large values of Ma would influence the temperature distribu- 
tion by convecting the warm liquid toward the center, whereas at lower values of Ma, similar to the 
values used in this study, no significant effect on temperature was observed. Therefore, J * is inversely 
proportional to Rs, i.e., J* °c l/Rs where Rs decreases with increasing AT. 

Another effect of combined flow is the influence of the dynamic pressure on the surface curva- 
ture. For the case of evaporation, the liquid is superheated with respect to the vapor and the thermocapil- 
lary stress causing the liquid to be driven out of the comer regions, thus causing a pressure decrease in 
this region. Conservation of mass causes the center portion of the surface to rise above that of the basic 
state. Figure 16 shows the divergence of the surface position from the basic state with increasing AT. 



Figure 16. Surface divergence from basic state for various AT. 


F. Combined Flow With Sidewall Temperature Gradient 

In order to evaluate the effect of a low thermal conductivity porous plug material, the sidewall 
boundary condition is adjusted to account for a temperature gradient through the porous plug. The mag- 
nitude of the gradient is determined by the Biot number, Bi, specified at the porous plugs upper solid 
boundary. The Biot number is defined as the ratio of convective heat transfer from the solid/vapor to 
conductive heat transfer from the liquid/solid. A low value of Bi specifies that the heat transfer is domi- 
nated by conduction, which is indicative of materials which have high thermal conductivities such as 
copper. Three values of Bi are used to evaluate this effect, Bi equal to 0.0, 0.15, and 0.5, and the corre- 
sponding sidewall temperature distributions are shown in figure 17. 

Using the above stated values of Bi, the resultant temperature and flow field is determined and 
the evaporative mass transfer calculated. Figure 1 8 shows a typical temperature and flow field for a 
porous plug with an imposed temperature gradient on the sidewall boundaries. 
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Figure 17. Sidewall boundary temperature for various Bi. 



Figure 18. Temperature and stream function distribution, Ma = 85.29, Bi =0.15. 

The effect of an imposed sidewall temperature gradient is most evident in the solution for the 
surface temperature distribution. Figure 1 9 shows the effect of Bi on the resultant surface temperature 
for A T = 1 . As Bi is increased, the heat transfer through the plug solid material is becoming dominated 
by convection on the vapor side, thus the surface remains cooler. This lowering of the surface tempera- 
ture results in a decrease in the evaporation across the interface as shown in tables 8 through 1 1 . 
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xl 

Figure 19. Surface temperature for various Bi, A T= 1.0. 

As shown in table 8, T sur f shows the same variation with AT as was observed in table 6, with 
increasing Bi causing a lowering of the curve or a decrease in T sur j. 

Table 8. T sur j versus Bi and AT. 


Biot No. 
(Bi) 

> 

ll 

o 

T 

1 surf 

AT = 1.25 

in 

ll 

h* 

< 

Basic 


- 

- 

0.00 



0.086 

0.15 

0.068 

| 

0.068 

0.50 

0.046 

0.046 

0.046 


Using the scaling relations given in equation ( 1 2) and the relevant properties given in table 1 , the 
dimensional values of T sur f are given in table 9. 

Table 9. T sur j versus Bi and AT (dimensional). 


Biot No. 
(Bi) 

O 

II 

< 

T SU rf (K) 
AT= 1.25 

in 

ll 

< 

Basic 

19.38 

- 

— 

0.00 

19.38 

19.15 

18.93 

0.15 

19.36 

19.13 

18.90 

0.50 

19.34 

19.10 

18.87 


Based on the surface temperature distribution given in tables 8 and 9, the resultant evaporation 
for the range of superheats considered are given in tables 10 and 1 1. 
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Table 10. 7* versus Bi and AT. 


Biot No. 

(Bi) 

> 

II 

© 

7* 

AT = 1.25 

AT = 1.5 

Basic 


- 

- 

0.00 

0.086 

0.107 

0.129 

0.15 

0.068 

0.086 

0.103 

0.50 

0.047 

0.061 

0.073 


Using the scaling relations given in equation (12) and the relevant properties given in table 1 , the 
dimensional values of 7 are given in table 1 1 . 

Table 11.7* versus Bi and AT (dimensional). 


Biot No. 
(Bi) 

> 

II 

o 

7 (kg/m 2 -s) 
AT = 1.25 

in 

II 

< 

Basic 

1.136 

- 

- 

0.00 

1.132 

1.410 

1.699 

0.15 

0.905 

1.132 

1.357 

0.50 

0.629 

0.804 

0.969 


Figure 20 shows the same linear trend of 7* versus AT across the pore as was shown in figure 15. 
The effect Bi is to lower the curve, while keeping the same general slope. 



Figure 20. 7* versus AT for various Bi. 

Figure 21 shows that the variation of evaporation flux, 7* as a function of Bi is nonlinear and 
decreases as Bi is increased. 
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Figure 21.7* versus Bi for various AT. 


G. Summary 

In this section, the relevant fluid properties and nondimensional groups for liquid hydrogen are 
presented. Superheats of 1 .0, 1 .25, and 1 .5 K are used in defining the nondimensional groups that are 
dependent on AT. To evaluate the effect of grid refinement on the solution, various grid sizes were 
examined. The parameters used to evaluate this effect are the comer surface temperature gradient, T , x , 

and the resultant evaporation flux across the pore. Grid sizes of 10x15, 15x20, 20x25, and 25x30 were 
used. It was found that as the grid is refined, i.e., the number of grid points is increased, both T, x and 7* 

are converging on their exact value. Although a grid size of 25x30 gave the smallest percent change 
from the pervious grid, a grid size of 20x25 is used in this study due to computer limitations. 

The steady state temperature distribution and flow field for the case of pure evaporation was 
solved and used as the basic state for comparison with the combined flow and sidewall thermal gradient 
boundary condition. The basic state was also compared to pure thermal conduction in the liquid. It was 
observed that the mean surface temperature, T sur j, decreased with evaporation due to the transfer of the 
latent heat across the interface. 

The steady-state temperature distribution and flow field for the case of combined evaporation/ 
thermocapillary flow was solved for various levels of liquid superheat. Pore temperature differentials, 
AT, of 1 .0, 1 .25, and 1 .5 K were used. It was shown that as the total evaporation across the pore, 7* , is 
directly proportional with AT. 

The steady-state temperature distribution and flow field for the case of combined evaporation/ 
thermocapillary flow with an imposed sidewall thermal gradient was solved using various values for the 
Biot number, Bi, i.e., the magnitude of convection from the pore solid boundary to the magnitude of 
solid conduction through the pore. Biot numbers of 0.0, 0.15, and 0.5 were used over the same range of 
AT used previously. As with the case of combined evaporation/thermocapillary flow with isothermal 
boundaries, 7* is directly proportional and linear with AT. The effect of Bi is to shift the curve down, 
i.e., as Bi increases the total evaporation across the pore decreases. This decrease in 7* is inversely 
proportional with Bi and nonlinear. 
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VI. DISCUSSION AND CONCLUSIONS 


A. Introduction 

The objective of this study was to investigate some of the physical parameters that influence the 
mass transport within a porous plug trapping device. In this study, the influence of parameters such as 
the level of superheat and the magnitude of an imposed thermal gradient along die sidewall was shown. 
The influence of other parameters such as Ma, Cr, Rs , Gr, and Vr were also indirectly investigated since 
all depend on the liquid superheat, i.e., AT across the pore. 

The effect of these parameters on the thermal distribution, flow field, and total evaporation is 
discussed. 


B. Discussion 

1. Basic State (Pure Evaporation) . The basic state as defined in this study is the temperature dis- 
tribution, flow field, and total evaporation in a pore having a diameter of 1 /im, a liquid superheat of 

1 K, and a vanishing Marangoni number. The flow field resulting from the basic state conditions is 
vertically oriented, i.e., the velocity vectors are directed from the bottom of the pore to the top and exit 
the pore through evaporation across the free surface. The only effect that drives the flow field in this 
case is the mass transport across the interface. Near the interface, the flow is directed toward the corners 
of the pore due to the high surface thermal gradient that occurs in the comers. Due to the high thermal 
gradient in the comer regions, the majority of the evaporation occurs in the corner areas with very little 
evaporation at the pore center. 

Liquid convection due to the mass transport across the interface provides surface cooling by 
driving the cooler liquid in the center of the pore to the warm comer regions and by latent heat transport 
at the liquid surface. 

2. Combined Evaporation/Thermocapillarv Flow . The flow field which results from the com- 
bined evaporation/thermocapillary driven flow consists of two counter-rotating cells near the liquid sur- 
face. The cells are formed by liquid convection along the surface caused by the variation of surface ten- 
sion. For this study, the surface tension is defined as being a decreasing linear function of temperature. 
Since the surface temperature decreases toward the center of the cavity, the liquid is convected toward 
the center. 

Thermocapillary convection will tend to raise the surface temperature by convecting warmer 
liquid from the pore corner regions to the cooler center, i.e., opposes the previously described evapora- 
tion effects where the cooler liquid in the center of the pore is convected toward the corners. As was 
shown in tables 6 and 7, as the AT across the pore increased, the liquid evaporation also increased while 
the mean surface temperature decreased. This would seem to indicate that die convection due to the 
evaporation and latent heat transport seemed to dominate the convection arising from thermocapillarity. 
Although, it was shown by Schmidt 8 that for very large values of Ma, i.e., Ma > 1,000, T sur j will tend to 
increase. 

3. Combined Flow With Sidewall Temperature Gradient . The flow field that results from 
combined flow with sidewall temperature gradient consists of two counter-rotating cells much 
like the flow field for isothermal sidewalls. The most significant difference relative to combined 
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evaporation/thermocapillary driven flow is the steady-state thermal distribution. The effect of a sidewall 
Bi is the formation of a thermal gradient throughout the body of the liquid. This is in contrast to the 
results discussed previously where the thermal gradient was restricted to the region about the free 
surface. Again, liquid is convected along the free surface and drives the warmer liquid from the corner 
regions to the cooler pore center. 

The most significant effect of Bi is the lowering of the overall surface temperature distribution, 
which is evident in both the actual temperature distribution across the pore and in T sur j. The surface 

temperature and temperature gradient in the corner regions decreases with increasing Bi. This reduction 
in temperature occurs since the heat conduction from the pore interior to the surface is reduced; thus, the 
sidewall temperature is lower than if the wall were isothermal, T l2 = 1.0. Since less heat is transported 

to the free surface, the overall evaporation across the interface is reduced accordingly. The variation of 
J* with AT is linear with the effect of Bi causing the curve to be lowered with increasing Bi. A secondary 
observation is that the variation of J* with Bi is nonlinear and decreasing. In this case, the effect of AT is 
to raise the curve with increasing AT. 


C. Conclusions 

The results of this study show that the performance of the porous plug type trapping device is 
influenced by external parameters, i.e., AT across the plug and the type of material that the plug is fabri- 
cated from. For plug type trapping devices, the two most important features are the maximum sustain- 
able internal pressure or bubble point and the influence of controllable and uncontrollable parameters on 
the evaporation flux across the plug. Controllable parameters are the pore diameter and length, number 
of pores, the plug material, and back pressure conditions. The uncontrollable parameters are the non- 
dimensional groupings that are dependent on the specific cryogen, e.g., liquid hydrogen. 

The analysis showed that the evaporation flux is strongly dependent on the temperature differen- 
tial across the pore and the pore surface temperature distribution. For Marangoni numbers representative 
of liquid hydrogen at the various superheated conditions, no significant effect on the surface temperature 
was observed. The evaporation flux increased primarily due to the reduction of Rs. It was also observed 
that the increase in Marangoni number caused the dynamic pressure in the corner regions to decrease 
due to liquid being driven toward the center. This causes the dynamic pressure in the center to increase, 
thus increasing the curvature of the surface and decreasing the effective bubble point of the pore. As 
shown in figure 17, the increase in surface curvature is small, thus resulting in a small decrease in bubble 
point. 


The effect of Biot number is most evident in the surface temperature distribution. As Bi 
increased, the surface temperature decreased, thus strongly influencing the evaporation across the pore. 
Although Bi caused a reduction in evaporation, it also reduced the strength of the surface thermal gradi- 
ent, which reduced the thermocapillary convection and, subsequently, the dynamic pressure within the 
pore. This lessened the effect of dynamic pressure on the effective bubble point of the pore. 

An important outcome of this study is the ability to numerically analyze the thermal distribution 
and flow field within the mircopores of the overall device. Previous empirical analysis on plug type 
trapping devices have measured the bubble point and evaporation rate for the specific plug tested on 
a macroscopic scale since no detailed analytical tool was available for performance predictions. As an 
example, assuming that a porous plug with a surface area of 4.9 cm 2 and a porosity of 30 percent was 
installed into the wall of a liquid hydrogen dewar. Assuming a AT = 1.0 K and accounting for 
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thermocapillary convection, the total evaporation flow rate across the plug would be approximately 

0.01665 gm/s. 


D. Recommendations 

1. Porous Plug Applications . Based on the results discussed, the design of the porous plug 
should be driven by the application. If maximum evaporation flux across the plug is desired, the use of a 
highly conductive material for the plug material is recommended. This reduces the side boundary Bi , 
thus transporting the maximum heat from the liquid to the free surface where evaporation is free to 
occur. On the other hand, if the plug is to be used as a phase separator for liquid acquisition, the use of a 
material that has a low thermal conductivity is recommended. This would serve to decrease the free sur- 
face thermal gradient, thus decreasing the thermocapillary convection and dynamic pressure and mini- 
mizing the decrease in bubble point. 

The use of uniform cylindrical pores for all applications is desirable due to the simple, pre- 
dictable, and repeatable geometry. This enables detailed microscopic analysis and allows the results to 
apply to similar plugs operated under identical conditions. 

2. Future Research . The analysis performed in this study was limited in several ways. The fol- 
lowing enhancements would add fidelity to the analysis. 

(a) The validation of these results through experimentation. Simple test cases can be designed 
which could serve to confirm the results presented in this study. Macroscopic tests can be used to verify 
the evaporation across the plug. With the number of pores known, the evaporation across each pore can 
be determined. To evaluate the How field distribution, a single, two-dimensional pore can be used with 
either liquid hydrogen or water with the relevant parameters adjusted by using chemical additives. The 
use of a microscope can be helptul in the flow field visualization. This will also serve to assess the accu- 
racy of this model. 

(b) The expansion of the analysis to three dimensions. This would more accurately model the 
actual physical geometry of the pore. 

(c) A closer and more extensive examination of the effect of pressure differential across the pore 
on the resulting temperature and flow fields. 

(d) Allowing the interface to translate within the pore depending on the applied pressure differ- 
ential and the energy balance across the pore. This would serve to address previous studies which sug- 
gest that the interface will retreat into the pore with increasing AT and evaporation, i.e., pore dry out. 
This would also require that the transient solution be solved. 


E. Summary 

In this section, the results for the case of pure evaporation, combined evaporation/thermo- 
capillary flow, and combined flow with a sidewall temperature gradient are discussed. The results seem 
to indicate that the material chosen for the fabrication of the plug should be application dependent. It 
was shown that a material with a high thermal conductivity, i.e., low Bi, serves to increase the free 
surface temperature, thus increasing the evaporation. Therefore, a plug fabricated out of a this material 
would be best used for an application where heat rejection or vapor generation is required. On the other 
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hand, if phase separation or liquid acquisition is required, a material with a high Bi should be used. This 
would serve to decrease the effect of the liquid dynamic pressure on the plug’s bubble point. 

In addition, several recommendations for the continuation of this study in terms of further analy- 
sis and experimentation are presented. 
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APPENDIX A 


FREE SURFACE SOLUTION METHOD 


A. Introduction 

In order to calculate the free surface position the nonlinear differential equation shown in section 
IV, equation (46), with the associated boundary conditions is solved using Newton’s method for non- 
linear systems. The methodology as it is applied to this problem is described in this section. 


5: 


B. Newton’s Method for Nonlinear Systems 

In order to calculate the free surface position, the coordinate system defined previously in figure 
x2, y 



xl=x- 0 xl-x= 1 


Figure A. 1. Free surface coordinate system. 

where x = xl, y = y (s) -x2, and y6) j s the surface coordinate referenced from the bottom of the cavity. 

Recalling the relations which defines the surface shape, namely, equation (46) and associated 
boundary conditions: 


--hTjK-'T’-*- 

(A.l) 

O 

II 

* 
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(A.2) 

CV,Jc)jc=i/2 “ 0 • 

(A.3) 
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and expressing equation (A.l) in finite difference form: 


>i+i ~ 2>y + y,-i 

h 2 


Boyi + II 
l-£> 



y«+i-y.--i 

2h 


»n3/2 


= 0 , 


(A.4) 


Newton’s method can be applied. 

The approach that is used to solve for the surface is to use Newton’s method for nonlinear 
systems to converge on a solution for the independent variable y/. Newton’s method employs the 
following iteration algorithm: 


• (a.5) 

where J t j is defined as the jacobian of F/, and is given by: 



\BF X 


BF,_ 

Byi 

dy 2 

dyj 

dF 2 

dF 2 


.By i 

By 2 . 


and F[ is defined as the root form of equation (A.4), namely: 


F,myM ' 2 » +^-1 


n3/2 


(A.l) 


where h is the separation distance between nodes, k is the iteration number and i,j are tensor index 
values. Expressing equation (A.l) in the general form: 


=y«+i -2y, + y M — )) 


(A.8) 


the Jacobian can be found. 

Equation (A.8) is differentiated with respect toy/ and the following relation for the Jacobian is 
obtained: 


/ _ dFj dy\ +l dyj dy,_ x 2 df 

,J Zyj ' 9yj By j % By j ' 


(A.9) 


Letting dy^dy^ = <5 1; , equation (A.9) can be written as. 
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(A. 10) 


Jij-fy+ij 28 jj + <5,_ ly - h ^ , 

where 5,- y is the Kronecker delta and is defined as having a valve of 1 when i=j and 0 when i*j. 


In order to determine df /dyj , f is expressed as: 



f(xi , y { , y'iiVi+iji-i )) = fix, , h(yi ),giy, )) , 

(A.l 1) 

where. 


h(y,) = y, , 

(A. 12) 

and 



(A. 13) 


2h 

Using equations (A.l 1), (A. 12), and (A. 13), dffdyj can be expressed as: 


ii 

S % 

(A. 14) 

_ df dy { df d ( y i+l - y,_ 1 3 
dh <hj dg dyj v 2 h ) 

(A. 15) 


(A. 16) 

and the expression for J,j can be written as: 



fj = 8 i+l j - 2 Sjj + 8i_ h j - h 2 


(A. 17) 

By inspection of equation (A. 17), J V) will be nonzero only in the range of i - 

-1<;</ + 1.. 

The nonzero values of A are given below: 



j = i ~ 1 : 

j..= i+ h % 
" 2 dy\ 

(A. 18) 

j = i ■ 

7=2 hdf 
U 2 dy, ’ 

(A. 19) 

j=i+\ : 

j.. = \---L. 

,J 2 dy\ ’ 

(A.20) 


where / is given by: 
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/ = 


Boy ; + n 
1 - Z) 



Differentiating / with respect to y ,• and yj , results in the following: 


and 



(A.22) 

(A.23) 


Now solving for the root equation over the entire domain, it is observed that at / = 1 and 
i = N + 1 , there are no values for j = i - 1 and 7 = 1 + 1, respectively, where N = the number of intervals 
in the domain. 

Considering the boundary conditions given in equations (A.2) and (A.3) for / = 1 and i = N + 1, 
respectively. Writing equation (A.7) at i = 2, and applying the boundary condition given in equation 
(A.2), the following is obtained: 


for i = 2: 


F 2 = y } - ly 2 + y, - h‘ 


Boy 2 + n 
1 - D 
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X3 — yi 

2h 
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(A.24) 


and at i =1: 

yi = 0 , 


(A.25) 


therefore, equation (A.24) becomes: 
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Now, considering the boundary condition at i = N + 1 : 
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(A. 27) 


and imposing the boundary condition at i = N + 1 : 

yN ± 2 zlE. = 0 


(A.28) 


or y N+2 = y N , equation (A.27) becomes: 
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(A. 29) 
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= 2 y N ~ 2y 
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Therefore, F, is defined over the interval of i = 2 to i = N+l. If one now shifts the indices over 
the interval of / = 1 and / = N, the following is obtained: 


at i = 1 : 
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at / = N: 


Fn - 2>’n-i ~ 2 yN ~ h 2 


Bny N + n \ 
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(A.32) 


Using the above definitions for the coefficients of the Jacobian, equations (A. 18), (A. 19), and 
(A.20), and the root equation given in equations (A. 30), (A.31), and (A.32), the independent variable y is 
solved for by iteration, using the algorithm given in equation (A.5). 


C. Summary 

In this appendix, Newton’s method for nonlinear equations that is used to solve for the free 
surface of the porous plug is discussed. Appendix B lists the numerical representation of this method. 
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APPENDIX B 


MODIFIED COMPUTER CODE 


A. Introduction 

In this section the modifications that were made to the reference model numerical method are 
shown. This section and reference 8 taken together represents the complete numerical model used in this 
study. 


B. Meniscus Subroutine 

SUBROUTINE MENISCUS(BO,CA,MA,PR,CR,RASPECT,FRHO, 

& VREC,NV AL,NXP,PULL,PV APJTYPE) 

£******************************************************************* 
C THIS PROGRAM WILL DEFINE A SURFACE INTERFACE BY SOLVING 
C A NON-LINEAR FINITE DIFFERENCE EQUATION BY USING NEWTONS 
C METHOD 

£********** *************************************************** ****** 
IMPLICIT DOUBLE PRECISION (A-H) 

IMPLICIT DOUBLE PRECISION (O-Z) 

CALCULATE ARRAY DIMENSIONS 

PARAMETER (MXE= 30, MYE=35) 

PARAMETER (MXN= 1+2*MXE) 

PARAMETER (MYN= 1+2*MYE) 

COMMON/SURFACE/ 

& XP(MXN),YP(MXN),YSURFO(MXN),PD(MXN),PB(MXN),TEMP(MXN), 
& SURFNN(MXN,2),SURFNT(MXN,2) 


DIMENSION ARRAYS AND DECLARE VARIABLES 

DIMENSION Y(MXN),TERM1(MXN),F(MXN),DEL(MXN),DY(MXN),R(MXN), 
& AJAC(MXN,MYN),TBAR(MXN),TERM2(MXN) 

REAL MAIN 

DEFINE FUNCTIONS 

FUNCl(X2,DYDX,BOND,CONl,CON2)=((BOND*X2+CONl)/(l.-CON2))* 

& (l.+DYDX**2.)**1.5 

FUNC2(DYDX, BOND ,CON2)=(BOND/(l.-CON2))*(l.+DYDX**2.)** 1.5 
FUNC3(X2,DYDX, BOND, CONI, CON2)=3*((BOND*X2+CONl)/(l.-CON2))* 

& SQRT(1.+DYDX**2.)*DYDX 
C 


PA GE B LANK NOT FILMED 


A 1 



DATA MAXITER/ 100000/ MAXIMUM NUMBER OF ITERATIONS 
DATA Yl/0./ [BOUNDARY VALUES 

DATA EANG/20./ [INITIAL ELEVATION ANGLE 

DATA STOL/O.OOO 1/ [SURFACE CONVERGENCE TOLERANCE 

C 

C INITIALIZATION OF VARIABLES 
C 

DELX=XP(2)-XP( 1 ) 

EANGR=EANG/57. 2957795 1 
C 

DO 10 N=1,NVAL 
R(N)=XP(N+l)/COS(EANGR) 

Y(N)=R(N)*SIN(EANGR) 

DO 20 M=1,NVAL 
AJAC(N,M)=0. 

20 CONTINUE 
10 CONTINUE 
C 

C DEFINE DELTA TEMPERATURE FOR THERMOCAPILLARITY TERM 
C 

DO 15 N=1,NVAL+1 

IF(ITYPE.EQ.0)THEN 

TBAR(N)=TEMP(N) 

ELSE 

TB AR(N)=TEMP(N)+ 1 . 

ENDIF 

15 CONTINUE 
C 

C LINEARLY INTERPOLATE ON PD 
C 

DO 16 N=4,NXP-3,2 

PD(N)=PD(N- 1 )+(((PD(N+ 1 )-PD(N- 1 ))/(XP(N+ 1 )-XP(N- 1 ))) 

& *(XP(N)-XP(N-1))) 

16 CONTINUE 
C 

C CALCULATE ROOT VECTOR 
C 

DELTA0=100. 

DO 30 ITER=1, MAXITER 
DO 40 N=1,NVAL 
IF(N.EQ. 1)THEN 

DY(N)=(Y(N+ 1 )-Y 1 )/2./DELX 
F(N)=Y(N+1)-2.*Y(N)-Y 1 
ELSEIF(N.EQ.NVAL)THEN 
DY(N)=0. 

F(N)=2.*(Y(N-1)-Y(N)) 

ELSE 
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DY(N)=(Y(N+ 1 )- Y(N- 1 ))/2./DELX 
F(N)=Y(N+1)-2.*Y(N)+Y(N-1) 



ENDIF 


C 

C CALCULATE SURFACE INFLUENCE DUE TO CA,PD,PULL,PVAP,VR,T,CR,TBAR 
C 

TERM1(N)=-CA*(PULL+PD(N+1)-PVAP)+VREC*TEMP(N+1)**2. 
TERM2(N)=CR*TBAR(N+ 1 ) 

C 

F(N)=F(N)-FUNC1(Y(N),DY(N),B0,TERM1(N),TERM2(N))*DELX**2. 

40 CONTINUE 

c 

C CALCULATE THE JACOBIAN 
C 

DO 50 N=I,NVAL 
IF(N.NE.l) AJAC(N,N- 1 )= 1 ,+DELX* 

& FUNC3(Y(N),DY(N),BO,TERM 1 (N),TERM2(N))/2. 

AJAC(N,N)=-2.-FUNC2(DY(N),BO,TERM2(N))*DELX**2. 

IF(N.NE.NVAL) AJAC(N, N+ 1 )= 1 .-DELX* 

& FUNC3(Y(N),DY(N),BO,TERM 1 (N),TERM2(N))/2. 

50 CONTINUE 
C 

C SOLVE FOR CHANGE IN Y-VALUE 
C 

CALL CROUT(AJAC,F,DEL,NVAL) 

C 

DELTA =0. 

DO 60 N=1,NVAL 
DELTA=DELT A+DEL(N ) 

Y(N)=Y(N)-DEL(N) 

60 CONTINUE 
C 

IF(ABS(DELTA).LE.STOL)GOTO 500 
DELTA0=DELTA 
C 

30 CONTINUE 

WRITE(8,4000) 

4000 FORMAT(2X, INTERFACE SOLUTION DID NOT CONVERGE',//) 

STOP 

500 CONTINUE 
C 

C DETERMINE X & Y FOR RIGHT SIDE OF INTERFACE 
C 

YP(1)=Y1 

DO 65 N=2,NVAL+1 
YP(N)=Y(N-1) 

65 CONTINUE 
K=0 

DO 70 N=NVAL+2,NXP 
K=K+1 

YP(N)=YP(N-2*K) 
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70 CONTINUE 

DO 75 N=1,NXP 
YP(N)=YP(N)+RASPECT 
75 CONTINUE 

DO 80 N=1,NVAL 

DYDX=(YP(N+1)-YP(N))/DELX 

ALFA=DATAN(DYDX) 

SURFNN(N, 1 )=- 1 . *DSIN( ALFA) 
SURFNN(N,2)=DCOS(ALFA) 
SURFNT(N,l)=DCOS(ALFA) 
SURFNT(N,2)=DSIN(ALFA) 

80 CONTINUE 

SURFNN(NV AL+ 1 , 1 )=0. 
SURFNN(NVAL+1,2)=1. 

SURFNT (NV AL+ 1 , 1 )= I . 

SURFNT (NV AL+ 1 ,2)=0. 

K=0 

DO 85 N=NV AL+2,NXP 
K=K+1 

SURFNN(N, 1 )=- 1 ,*SURFNN(N-2*K, 1 ) 
SURFNN(N,2)=SURFNN(N-2*K,2) 
SURFNT(N,1)=SURFNT(N-2*K,1) 
SURFNT(N,2)=- 1 ,*SURFNT(N-2*K,2) 
85 CONTINUE 

RETURN 

END 


C. Crout Subroutine 


SUBROUTINE CROUT(A,B,X,N) 

(2 * * * ** * * * * * * * **************** * * * ** * ************ **** * * * ***** *** * **** * 
C SOLVES FOR THE X MATRIX USING A CROUT REDUCTION OF A NXN 
C TRIDIAGONAL MATRIX, AX=B 

£************* ****** ******************************************* ****** 
IMPLICIT DOUBLE PRECISION (A-H) 

IMPLICIT DOUBLE PRECISION (O-Z) 

C 

PARAMETER (MXE= 30, MYE= 35) 

PARAMETER (MXN= 1+2*MXE) 

PARAMETER (MYN= 1+2*MYE) 

C 

DIMENSION A(MXN,MYN),B(MXN),X(MXN),C(MXN,MYN), 

& U(MXN,MYN),Z(MXN) 

C 

NM1=N-1 


FIRST ROW 
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c 

C(1,1)=A(U) 

U(1,2)=A(1,2)/C(1,1) 

C 

C ROWS BETWEEN 1 AND N 
C 

DO 10 1=2, NM1 
C(I,I-1)=A(I,I-1) 

C(I,I)= A(I,I)-C(I,I- 1 )*U(I- 1 ,1) 
U(I,I+1)=A(I,I+1)/C(I,I) 

10 CONTINUE 
C 

C LAST ROW 
C 

C(N,N- 1 )=A(N,N- 1 ) 
C(N,N)=A(N,N)-C(N,N- 1 )*U(N- 1 ,N) 
C 

C SOLVE FOR LZ=B 
C 

Z(1)=B( 1 )/C( 1,1) 

DO 20 1=2, N 

Z(I)=(B(I)-C(I,I- 1 )*Z(I- 1 ))/C(I,I) 

20 CONTINUE 
C 

C SOLVE FOR UX=Z 
C 

X(N)=Z(N) 

DO 30 1=1, NM1 
K=N-I 

X(K)=Z(K)-U(K,K+1)*X(K+ 1 ) 

30 CONTINUE 
C 

RETURN 

END 


D. Sidewall Boundary Condition Modification 

C 

C SET LEFT AND RIGHT BOUNDARY CONDITIONS 
C 

IF(ITBC.EQ.0)THEN 

NL=1 

NR=NXP 

DO 157 N=1,NYP 

TEMP(NL)=TSF( 1 ) 

TEMP(NR)=TSF( 1 ) 

NL=NL+NXP 
NR=NR+NXP 
157 CONTINUE 
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TEMP(NGNODE+5)=TSF( 1 ) 

TEMP(NGNODE+ 1 2+NMID)=TSF( 1 ) 

TEMP(NGNODE+ 1 7+NMID)=TSF( 1 ) 

TEMP(NGNODE+24+2*NMID)=TSF(l) 

ELSEIF(ITBC.EQ. 1 )THEN 

NL=1 

NR=NXP 

DO 158 N=1,NYP 

TEMP(NL)=TSF(2)-NU*(TSF(2)-TVAP)/(1.+NU*RASPECT)*REAL(XG(NL,2)) 

TEMP(NR)=TEMP(NL) 

NL=NL+NXP 

NR=NR+NXP 

CONTINUE 

TEMP(NGNODE+5)=TSF(2)-NU*(TSF(2)-TVAP)/(l.+NU*RASPECT) 

*REAL(XG(NGNODE+5,2)) 

TEMP(NGNODE+ 1 2+NMID)=TEMP(NGNODE+5) 

TEMP(NGNODE+17+NMID)=TSF(2)-NU*(TSF(2)-TVAP)/(l.+NU*RASPECT) 
*REAL(XG(NGNODE+ 1 7+NMID.2)) 

TEMP(NGNODE+24+2*NMID)=TEMP(NGNODE+24+2*NMID) 

ENDIF 


E. Lower Boundary Condition Modification 


UM AX=(3./2. )*QFLOW 

ACAV=0.5 

DO 69 NBOT=l,NXP 

VELBOT(NBOT)=UMAX*(l.-((REAL(XSURF(NBOT))-ACAV)**2./ACAV**2.)) 

CONTINUE 
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